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Abstract 

We study the non-equilibrium dynamics of a pair of qubits made of two-level atoms separated 
in space with distance r and interacting with one common electromagnetic field but not directly 
with each other. Our calculation makes a weak coupling assumption but no Born or Markov 
approximation. We derived a non-Mar kovian master equation for the evolution of the reduced 
density matrix of the two-qubit system after integrating out the electromagnetic field modes. It 
contains a Markovian part with a Lindblad type operator and a nonMarkovian contribution, the 
physics of which is the main focus of this study. We use the concurrence function as a measure 
of quantum entanglement between the two qubits. Two classes of states are studied in detail: 
Class A is a one parameter family of states which are the superposition of the highest energy 



\I) = |11) and lowest energy \0) = |00) states, viz, \A) = ^/p\I) + -p)\0), with < p < 1; 
and Class B states \B) are linear combinations of the symmetric |+) = "^(|01) + 1 10) ) and the 
antisymmetric |— ) = — ^ ( 1 01} — |10)) Bell states. We obtain similar behavior for the Bell states 
as in earlier results derived by using the Born-Markov approximation [40] on the same model. 
However, in the Class \A) states the behavior is qualitatively different: under the non-Markovian 
evolution we do not see sudden death of quantum entanglement and subsequent revivals, except 
when the qubits are sufficiently far apart. (The existence of sudden death was first reported for 



two qubits in two disjoint cavity electromagnetic fields [38j, and the dark period and revival were 
found from calculations using the Born-Markov approximation For an initial Bell state, 

our findings based on non-Markovian dynamics agree with those obtained under the Born-Markov 
approximation. We provide explanations for such differences of behavior both between these two 
classes of states and between the predictions from the Markov and non-Markovian dynamics. We 
also study the decoherence of this two-qubit system and find that the decoherence rate in the case 
of one qubit initially in an excited state does not change significantly with the qubits separation 
whereas it does for the case when one qubit is initially in the ground state. Furthermore, when the 
two qubits are close together, the coherence of the whole system is preserved longer than it does 
in the single qubit case or when the two qubits are far apart. 
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I. INTRODUCTION 



Investigation of quantum entanglement is both of practical and theoretical significance: 
It is viewed as a basic resource for quantum information processing (QIP) ll| and it is a basic 
issue in understanding the nature of nonlocality in quantum mechanics 3], 4|. However, 
even its very definition and accurate characterization are by no means easy, especially for 
multi-partite states (see, e.g., 0, [g, Q, s|, [9], 



10 



111].) Nonetheless there are use 



proposed for the separability of a bipartite state, pure and mixed (e.g., 



12, 



13, 



14 



ul criteria 



15 



181]). theorems proven (e.g, 



s introduced (e.g. 



21 



IT 



23|. 



l20j), and new mathematical too 
which add to advances in the last decade of this new endeavor 

Realistic quantum systems designed for QIP cannot avoid interactions with their envi- 
ronments, which can alter their quantum coherence and entanglement properties. Thus 
quantum decoherence and disentanglement are two essential obstacles to overcome for the 
design of quantum computers and the implementation of QIP. Environment-induced deco- 
herence in the context of QIP has been under study for over a decade [24| and studies of 
environment-induced disentanglement has seen a rapid increase in recent years. There are 
now experimental proposals to measure finite time disentanglement induced by a reservoir 



251 ] . The relation between decoherence and disentanglement is an interesting one because 
both are attributable to the decay of quantum interference in the system upon the interaction 
with an environment. (See, e.g., 26. 27. 28. H. I30I]) 



In addition to the mathematical investigations mentioned above which could provide 
rather general characterizations of quantum entanglement, detailed studies of physical mod- 
els targeting actual designs of quantum computer components can add precious insight into 
its behavior in concrete settings. Two classes of models relevant to condensed matter and 
atomic-optical QIP schemes are of particular interest to us. The first class consists of the 
quantum Brownian motion model (QBM) and the spin-boson model (SBM). Quantum deco- 
herence has been studied in detail in both models, and results on quantum disentanglement 



are also appearing (See [26 



27 



23, 



291 ] for QBM under high temperature, negligible dissi- 



pation, and 3l| for an attempt towards the full non-Markovian regimes.) The second class 
of models describes atoms or ions in a cavity with a quantum electromagnetic field at zero 
or finite temperature. The model consists of two two level atoms (2LA) in an electromag- 
netic field (EMF). For a primary source on this topic, read, e.g.. 321 ] . For a more recent 
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description of its dynamics under the Born-Markov approximation, see the review of 



33] 



A. Two-atom entanglement via an electromagnetic field 



Quantum decoherence and entanglement between one 2LA and an EMF has been treated 
by us and many other authors earlier 



34 



351 ] and by Cummings and Hu recently towards the 



strong coupling regime [36], which provide insight in how the atom- field interaction affects 
their entanglement. There is recent report of exact solutions found for a 2LA in an EMF 

n 

using the underlying pseudo-supersymmetry of the system [37[ • In the 2 atom-EMF model, 
the two atoms can be assumed to interact only with its own cavity EMF, or with a common 
EMF, and they can also interact with each other. The noninteracting case of separate fields 



was first studied by Yu and Eberly 



38 



391 ] where 'sudden death' of quantum entanglement 



was si ghte d. The noninteracting case of a common field was studied recently by Ficek and 
Tanas 40], An et al \&^. Quantum decoherence of N-qubits in an electromagnetic field was 
studied by Palma, Suominen and Ekert 42] . For entanglement of ions in cavities, see, e.g., 



43]. 



For the purpose of quantum information processing, we have emphasized earlier in our 
studies of quantum decoherence that it is absolutely essential to keep track fully of the 
mutual influence of, or the interplay between, the system and the environment. If one 
chooses to focus only on the system dynamics, one needs to take into consideration the 
back-action of the environment, and vice versa. 
In our prior work 



34 



35], we used the influence functional formalism with a Grassman- 



nian algebra for the qubits (system) and a coherent state path integral representation for the 
EMF (environment). Here, we employ a more standard operator method through pertur- 
bation theory, because the assumption of an initial vacuum state for the EMF allows a full 
resummation of the perturbative series, thus leading to an exact and closed expression for 
the evolution of the reduced density matrix of the two qubits. This approach incorporates 
the back-action of the environment on the system self-consistently and in so doing generates 
non-Markovian dynamics. We shall see that these features make a fundamental difference 
in the depiction of evolution of quantum entanglement in the qubit dynamics. 
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B. The importance of including back-action self-consistently 



Since quantum entanglement is a more delicate quantity to track down than decoherence, 
an accurate description is even more crucial. For this, one needs to pay extra attention to 
back-actions. For example, in the case of two 2LA (system) in a cavity EMF (environment), 
the two parties are equally important. This means that we should include both the back- 
action of the field on the atoms, and the atoms on the field. In a more complete treatment 
as attempted here, we obtain results qualitatively different from earlier treatments where 



331 ] . Some special effects like 



the back-action is not fully included or properly treated 
'sudden death' [38j can in this broader context be seen as consequences only of rather special 
arrangements: Each atom interacting with its own EMF precludes the fields from talking 
to each other and in turn cuts off the atoms' natural inclination (by the dictum of quantum 
mechanics) to be entangled. In effect, this is only a limiting case of the full dynamics we 
explored here for the two-qubit entanglement via the EMF. This limit corresponds to the 
qubits being separated by distances much larger than the correlation length characterizing 
the total system. For a wide range of spatial separations within the correlation length, 
entanglement is robust: Our results for the full atom-field dynamics reveal that there is no 
sudden death. 



C. Non-Markovian dynamics from back-action 



It is common knowledge in nonequilibrium statistical mechanics |44| that for two inter- 
acting subsystems the two ordinary differential equations governing each subsystem can be 
written as an integro-differential equation governing one such subsystem, thus rendering its 
dynamics non-Markovian, with the memory of the other subsystem's dynamics registered 
in the nonlocal kernels (which are responsible for the appearance of dissipation and noise 
should the other subsystem possess a much greater number of degrees of freedom and are 
coarse-grained in some way). Thus inclusion of back-action self-consistently in general en- 
genders non-Markovian dynamics. Invoking the Markov approximation as is usually done 
may lead to the loss of valuable information, especially for quantum systems. These assump- 
tions need to be scrutinized carefully with due consideration of the different time scales in 
the system and the specific physical quantities that are of interest in the study. 
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For monitoring the evolution of quantum entanglement which is usually a more delicate 
process than coherence, if one lacks detailed knowledge of how the different important pro- 
cesses interplay, our experience is that it is safer not to put in any ad hoc assumption at the 
beginning (e.g., Markovian approximation, high temperature, white noise) but to start with 
a first principles treatment of the dynamics (which is likely non-Markovian) involving all 
subsystems concerned and respecting full self-consistency. This is because entanglement can 
be artificially and unknowingly curtailed or removed in these ad hoc assumptions. What is 
described here is not a procedural, but a substantive issue, if one seeks to coherently follow 
or manipulate any quantum system, as in QIP, because doing it otherwise can generate 
quantitatively inaccurate or even qualitatively wrong results. 

Thus the inclusion of backreaction (which depends on the type of coupling and the features 
of the environment) usually leads to nonMarkovian dynamics {|3lj]. Also, under extreme 
conditions such as imposing infinite cutoff frequency and at high temperatures, the dynamics 
of, say, a quantum harmonic oscillator bilinearly coupled to an Ohmic bath can become 
Markovian 45]. Other factors leading to or effecting nonMarkovian behavior include the 
choice of special initial conditions. For example, the factorizable initial condition introduces 
a fiducial or special choice of time into the dynamics which destroys time-homogeneity. 

A word about terminology might be helpful here: One usually refers to Markovian dynam- 
ics as that governed by a master equation with constant-in-time coefficients, i.e., described 
by a Linblad operator, and non-Markovian for all other types of dynamics. A more restricted 
condition limits the definition of non-Markovian dynamics to cases with non-trivial (nonlocal 
in time) integral kernels appearing in the master equation. This more stringent definition 
would refer to dynamics (depicted by master equations containing coefficients which are) 
both time-homogeneous and non-time-homogeneous as Markovian. We use the first and 
more common convention of terminology, in which the master equation (I85I) which is local 
in time but non-time-homogeneous would be nonMarkovian. The Markovian regime emerges 
in the limit when the two qubits are far separated. This feature is similar to the HPZ master 
equation 46J for quantum Brownian motion, where Markovian (time-homogeneous) dynam- 
ics appears only in specific limiting conditions (high temperature and ohmic distribution of 
environmental modes, as alluded to above). 

Our present study of the two qubit (2qb)- EMF system is also aimed at addressing a 
common folklore, namely, that in quantum optics one does not need to worry about non- 
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Markovian effects. We will see that there is memory effect residing in the off diagonal 
components of the reduced density matrix for the 2 qubit system which comes from virtual 
photon exchange processes mediated by the field and which depends on the qubit separation. 
Perhaps the simplest yet strongest reason for the need to take nonMarkovian effects seriously 
is that results from the Markovian approximation are incomplete and lead to qualitatively 
wrong predictions. 



D. Relation to prior work and organization of this paper 

In this paper we study the non-Markovian dynamics of a pair of qubits (2LA) separated 
in space by distance r interacting with one common electromagnetic field (EMF) through a 
Jaynes-Cummings type interaction Hamiltonian. We use the concurrence function as a mea- 
sure of quantum entanglement between the two qubits. The same model was studied before 



in detail 



3y Ficek and Tanas [33[ using the Born-Markov approximation. In a more recent 
paper 40| they show the existence of dark periods and revival of quantum entanglement in 
contrast to the findings of Yu and Eberly which assumes two qubits in disjoint EMFs. 

Our calculation makes a weak coupling assumption and ignores the two- and higher- 
photon- exchanges between the qubits, but it makes no Born or Markov approximation. We 
derive a non-Markovian master equation, which differs from the usual one of the Lindblad 
type: it contains extra terms that correspond to off-diagonal elements of the density matrix 
propagator. We concentrate on two classes of states, superpositions of highest and lowest 
energy states and the usual antisymmetric |— ) and symmetric |+) Bell states [3j and observe 
very different behavior. These are described in detail in the Discussions section. The 



difference between our results and that of Ref. [40j highlights the effect of non-Markovian 

(with memory) evolution of quantum entanglement. In short, we find similar behavior in 

the Class B (Bell) states but qualitative different behavior in the evolution of Class A states. 
□ 

Ref [40] found that their evolution leads generically to sudden death of entanglement and a 
subsequent revival. In our more complete treatment of the atom-field dynamics we indeed 
see the former effect present for large values of the inter-qubit distances. However, sudden 
death is absent for short distances, while there is no regime in which a revival of entanglement 
can take place. This calls for caution. 



Another set of papers close to our work reported here is that of 47] who considered two 
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2LA in an infinite temperature field bath. When the atoms are separated at large distance 
the authors assume that they are located inside two independent baths. ( The severance of 
the field is subject to the same criticism above: A small but finite quantum entanglement 
cannot be equated to zero because the small amount can later grow.) For these conditions 
and under the Markovian approximation, the time evolution of the two-atom system is given 
by the ergodic dynamical semi-group. They ignore without justification the effect of distance 
on the interaction between the qubits. A paper of interest not directly related to the present 
model but which does show the dependence of the disentanglement rate on distance, like 
ours reported here, is that by Roszak and Machinikowski 48J. They consider a system of 
excitons with different coupling, and with a very different infrared behavior of the bath 
modes. The latter seems not to be relevant to the two atoms' case here. 

This paper is organized as follows: Section 2 contains the main derivation. We write down 
the Hamiltonian for two 2-level atoms (2LA) interacting with a common electromagnetic 
field (EMF) at zero temperature, and we compute the relevant matrix elements for the 
propagator of the total system by resummation of the perturbative series (Appendix A). We 
then determine the evolution of the reduced density matrix of the atoms, which is expressed 
in terms of seven functions of time. We compute these functions using an approximation that 
amounts to keeping the contribution of the lowest loop order for the exchange of photons 
between the qubits. In Section 3 we examine the evolution of the reduced density matrix for 
two classes of initial states. We then describe the time evolution of quantum entanglement 
with spatial separation dependence in these states via the concurrence plotted for some 
representative cases. In Section 4, we study the decoherence of this system when the two 
qubits are initially disentangled. We consider two cases that correspond to one of the qubits 
being initially in a vacuum state and in an excited state. We compare these results with the 
single qubit cases and highlight the lessening of decoherence due to the presence of other 
qubit(s). In Section 5 we discuss and compare our results on disentanglement with the work 
of Yu and Eberly for two 2LA in separate EMF baths, and with the work of Fizek and Tanas 
on two 2LA in a common EMF bath under the Born-Markov approximation. We identify 
the point of departure of quantum dynamics under the Markovian approximation from the 
full non-Markovian dynamics and thereby demonstrate the limitations of the Born-Markov 
approximation. Finally, we discuss the domain of validity of the rotating wave approximation 
in describing these systems. In Appendix C we sketch an alternative derivation through the 



S 



Feynman- Vernon influence functional technique, in which Grassmann variables are employed 
for the study of the atomic degrees of freedom. 

II. TWO-ATOMS INTERACTING VIA A COMMON ELECTROMAGNETIC 
FIELD 

A. The Hamiltonian 

We consider two 2-level atoms (2LA), acting as two qubits, labeled 1 and 2, and an 
electromagnetic field described by the free Hamiltonian H 

H = hJ2 ^kfejA + nwoSPs™ + huoS^ s {2) (1) 

k 

th 

where is the frequency of the k l electromagnetic field mode and u Q the atomic frequency 
between the two levels of the atom, assumed to be the same for the two atoms. The 
electromagnetic field creation (annihilation) operator is (6 k ), while (S^) are the 
spin raising (lowering) operators for the rfi 1 atom. We will define the pointing vector from 
1 to 2 as r = r2 — ri and we will assume without loss of generality that ri + Y2 = 0. 

The two 2LAs do not interact with each other directly but only through the common 
electromagnetic field via the interaction Hamiltonian 

Hi = hJ29k (bl(e- ik - r/2 S { y + e^/ 2 ^) + b k (e^ 2 S? + e^ 2 ^)) , (2) 

k 

where g k = ^==, A being the coupling constant. We have assumed that the dipole moments 
of the atoms are parallel. The total Hamiltonian of the atom-field system is 

H = H + Ej. (3) 

B. Perturbative expansion and resummation 

We assume that at t — the state of the combined system of atoms+field is factorized 
and that the initial state of the EMF is the vacuum \0). For this reason we need to identify 
the action of the evolution operator e~ tHt on vectors of the form \0) <E> 1^)) where is a 
vector on the Hilbert space of the two 2LA's. 
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For this purpose, we use the resolvent expansion of the Hamiltonian 

dEe~ iEt 



-iHt 



E — H + irj 



(4) 



and we expand 



(E - H)- 1 — (E — Ho)' 1 + (E - H Q )- l Hj{E - H 



.-1 



{E - H y 1 H I (E - H y 1 H I (E - Ho)" 1 + ... 



(5) 



Of relevance for the computation of the reduced density matrix of the two qubits are 
matrix elements of the form (z; i',j'\(E — H)~ l \0] i,j), where z refers to a coherent state of 
the EM field and i,j = 0, 1, the value i = corresponding to the ground state of a single 
qubit and % — 1 to the excited state. We compute the matrix elements above through the 
perturbation expansion (J5]). It turns out that we can effect a resummation of the perturbative 
series and thus obtain an exact expression for the matrix elements-see Appendix A for details 
of the resummation. 

The non-vanishing matrix elements are the following 



kr 



~ (E — uj — a(E) - f3{E,r)e lkr ){E - uj k ) 



(6) 



>z;0,l\(E-H)- L \O;0,l) 



1 



1 



+ 



(z;l,0\(E-H)~ L \O;0A) 



E- u -a{E) - (3{E,r) 
1 

E-n-a(E)+P(E,r) 

1 T 1 

2 [E -uj -a{E) - f3(E,r 

1 



(7) 



E-iu -a(E)+ f3(E,r) 
[z;0,0\(E - Hy^O; 0,0} = E' 1 



{z- 1,1\{E - HY^O- 1,1) 

(z^O^E-HY^O^l) 

{z- 0,1\{E - HY^O- 1,1) 
{z-AME-HY'lO;!,!) 



1 



E - 2uo - 2a(E - uo ) - f(E,r) 



E 



k - E-2tu - 2a(E - uo ) - f(E, r 



E 



' e 2 



kk ^ (E - 2uj )(E -uj - ujy 
In the equations above the functions a(E), (3(E,r) are 

= E 



V 



e 2 



(8) 

(9) 
(10) 

(11) 

:i - L)£{12) 



E - a> k 



(13) 
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k h, — 

The definitions of the kernel if kk / and of the function / involve complicated expressions. 
However, the term involving i^kk' does not contribute to the evolution of the reduced density 
matrix, while the function f(E,r) is of order A 4 and it can be ignored in the approximation 
we effect in Sec. II. E. Thus for the purpose of this investigation, the explicit definitions of 
H and / are not needed, and hence not given here. 

Finally, the matrix L is defined as 



L : = 



'E 



e e 



(15) 



where E and are matrices on the space of momenta 

Hkk ' = E-ll-^ { a(E ~ " k)4k ' + 9k9k ' ( B^ + E-J-j ) (16) 

e - = e^t-j: ( i3(B - + ^ + e-^-J ) • (17) 

and the overbar denotes complex conjugation. 



C. The matrix elements of the propagator 



The next step is to Fourier transform the matrix elements of the resolvent in order to 
obtain the matrix elements of the evolution operator. Explicitly, 



(z;0,0|e-^|O;0,l) = £ e^ 2 ^(t) 



(z;0,l|e-^|O;0,l> = J 



dEe 



iEt 



1 



E-u -a(E)-(3(E,r) 



+ 



1 



E-LU -a(E)+[3(E,r)_ 



-■■ v + (t) 



(z;l,l|e-^|O;0,l) = J 



dEe 



iEt 



1 



E-u -a(E)-P(E,r) 



E-w -a(E)+P(E,r) 
(z;0,0\e- iAt \O;0,0) = 1 

(z;l,l|e-^|0;l,l) = J — 



■ v-(t) 



dEe 



-iEt 



2cu - 2a(E - cu ) - f(E,r) 



=: u(t) 



(18) 



(19) 



(20) 

(21) 
(22) 
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7 ^ E - 2u) Q - 2a(E - u ) - f{E,r) 



( 



'z;0,l\e- iHt \O;l,l) 



E- 



\ 



{z;l,0\ e - mt \O;l,l) ) V ~ ^ e^^(t) 
where we defined the functions Sk(i), ^k(^), ^k(^) as 



(23) 
(24) 



(*) = / 



(£ _ Wo _ a (E) - /3(E, r)e^)(E - uj k ) 

( 9k' 



E— U) — LU k l 

. 9k' 

\ E—0J o — Uly 



(25) 
(26) 



D. The reduced density matrix 

We next compute the elements of the reduced density matrix for the qubit system by 
integrating out the EM field degrees of freedom 

p&(*)= E p;g(o)/H[^^ (27) 

where [dz] is the standard Gaussian integration measure for the coherent states of the EM 
field. 

Substituting Eqs. (I18H24I) into (1271) we obtain through a tedious but straightforward 
calculation the elements of the reduced density matrix 



p^t) = P m\u\ 2 (t) (28) 

Poi(t) = p£(0)u(t)v* + (t) + p{l(0)u(t)v*_(t) (29) 

pll(t) = ^(o)«(t)<(t) + p£i(o)«(*)«-(*) (so) 

PooW = P>(t) (31) 

PooW = Poo(OKW + pMv-(t) + Poi(0)/i!(t) + d£(0)/ia(t) (32) 

pooW = pi>w(t) + f${o)v.(t) + pimm + pimm (33) 

p°o\(t) = P^(o)|^| 2 (t) + p?J(oK(tK(t) + P i°(o)|^| 2 (t) 

+ p&{0)v-{t)i%(t) + fftWKtf) (34) 
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+ f%(0)v-(t)vl(t) + f$(0)K*(t) (35) 
P§8(*) = 1 - Pn(t) - P°ol(t) - P \°o(t) (36) 



where 





= X)«k^( f ) s k(*) 
k 


/i 2 (t) 


k 




= EkTO 

k 




k 



(37) 

-ik-r ( 3g ) 

(39) 

k ' r , (40) 



and the functions w(t),v±(t) were defined in Eqs. f l2~2"j) . (ITU1) and (J2Uj) . 



E. Explicit forms for the evolution functions 

Eqs. (I28H36P provide an exact expression for the evolution of the reduced density matrix 
for the system of two qubits interacting with the EM field in the vacuum state. The evolution 
is determined by seven functions of time u, v±, Ki t 2, /xi,2, for which we have provided the full 
definitions. To study the details of the qubits' evolution we must obtain explicit forms for the 
functions above. For analytic expressions, we recourse to an approximation: Assuming weak 
coupling (A 2 << 1), we ignore the contribution of all processes that involve the exchange of 
two or more photons between the two qubits. 



1. The functions u, v± 



With the approximation above, the contribution of the function / drops out from the 
definition of u. Thus we obtain 

dEe~ iEt 



u(t) 



v± 



E — 2lu q — 2a(E - lu ) 



dEe 



-iEt ■ 



± 



1 



(41) 
(42) 



_E - lu - a(E) - P(E, r) E — lu q — a(E) + @(E, r) 

We evaluate these expressions using an additional approximation. In performing the 
Fourier transform, we only keep the contribution of the poles in the integral and ignore 
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that coming from a branch- cut 
expression of a(E)-see Ref. |34| 



hat appears due to the presence of a logarithm in the exact 
for details. We then obtain 



u(t) 
v±(t) 



-2iu) o t-2V t 
„ -iuiot—Fot 



(43) 
(44) 



In the equations above, we effected a renormalization of the frequency uj q by a constant 
divergent term- see [34] . The parameters r ,r r and a{r) are defined as 



T := —Ima{ui ) 
-a(r) + iT r := /3(u ,r), 



(45) 
(46) 



and they read explicitly 



r 
r,. 



A 2 u; 

A 2 sin u; r 



2nr 



A 2 



r 7T 



2vr 2 r 

+ sintt;o^[log(e 7 ct;o^ 



cosu r[- - Si(u r)\ 



w ° r , 1 — cosz. 


2 



(47) 
(48) 

(49) 



where Si is the sine-integral function. 

The term a(r) is a frequency shift caused by the vacuum fluctuations. It breaks the 
degeneracy of the two-qubit system and generates an effective dipole coupling between the 
qubits. At the limit r — >• 0, this term becomes infinite. One should recall however that the 
physical range of r is always larger than clb, the Bohr radius of the atoms. As r — > oo, 
a(r) -> 0. 

The constant T corresponds to the rate of emission from individual qubits. It coincides 
with the rate of emission obtained from the consideration of a single qubit interacting with 
the electromagnetic field. The function T r is specific to the two-qubit system. It arises from 
Feynman diagrams that involve an exchange of photons between the qubits. Heuristically, 
it expresses the number of virtual photons per unit time exchanged between the qubits. As 
such, T" 1 is the characteristic time-scale for the exchange of information between the qubits. 
As r — > 0, T r — > T and as r —>■ oo, T r — > 0. Note that the ratio r/r 
than unity, is of the order of unity as long as r is not much larger than oo^ 1 . 



smuJ ° r while smaller 
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It is interesting to note that T r = for r = mttUQ 1 , where n an integer. This is a resonant 
behaviour, similar to that of a classical oscillating dipole when r = n\/2, where A is the 
oscillation wavelength. 

2. The functions Ki,2(i) 



We first compute the functions z/ k , j/ k of Eq. (|26|) keeping terms up to second loop order 

/ .. U\ \ - m -iEt ( 




dEe 



E 



9k! 



E — 2u k , E — u — low 

where H and 6 are given by Eqs. (fT6j) and (JTTj). 

The summation over k' yields within an order of A 5 



<5kk' 

<5kk' 




(50) 




dEe 



iEt 



9k 



E — 2u n E 



UJ n 



a(E - cu k ) + (3(E - u k ) 



E -UJ n 



1 + 2 



a(E - u ) 



E-2u n 



x 



l + Ek' 



g?,(l +e i(k-k')-r) 



I — 1 — S . . 

1 , V g2 ;(1+e -i(k-k'). r) 

\ 1 + z^k' (E-w k -w k ,)(£;-w -w k ,) 



(51) 



The terms in brackets in the first line of the equation above can be absorbed in the 
leading-order denominators-see Eq. (152j) . The term in the second line, however, only gives 
rise to a (time- independent) multiplicative term of the form 1 + 0(A 2 ). Hence, if we keep 
the leading order terms in the expression of u^, we may ignore this term. Then z/k = f k , 
and within an error of order A 5 

dEe~ iEt 



Vk(t) =9k[ 



[E - 2u -2a(E- w )] [E - oo - cu k - a(E - cu k ) - (3(E - cu k , r))] 
Using the same approximation as in Sec. II.E.l for the Fourier transform, we obtain 



(52) 



Vk(t) = gk- 



-int-r t 



,—ioJot—Tot 



— e 



-iuj^t—iat—T 



uj — uj-k — a — ir + iT r 
We then substitute the expression above for into Eqs. (1391) and (|4"01) to get 
A 2 



(53) 



Ki(t) 



2vr 2 



2i'"' / °° uau e_2r °* + e ~ 2Vrt ~ 2e- (r ° +r '' )t cos[(cj — k — a)t] 

(3 I KiCLKi ■ 







(k-uj + af + (r - r r 



n 2 {t) 



A 2 



27r 2 r 



,oo e -2r t + e -avt _ 2 -(r +rv)t cog^ _ & _ a ui 
~ 2Tot 1 dk ^ sin Ax 



o 



(k-cu + af + (r - r r ) 2 



(54) 
(55) 
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For u t » 1, it is a reasonable approximation to substitute the Lorentzian with a delta 
function. Hence, 



K 1 (t) = r Q K(t) (56) 

K 2 (t) = T r K(t) (57) 



where 

1 



K{t) ~ — ^—e- 2r °*(e- rot - e-^) 2 - (58) 
r — r r 



3. The functions (j,^ 



We first compute the functions s^{t) of Eq. fl25|) 



g— iu) t— Tot— (r r +icr)e lk ' r t g— iui-^t 

= ^—7 ^TliTT- ( 59 ) 

Substituting into Eqs. fl37|) and fl38l) we obtain 

\2 ~i ~ ~—iuJot—Fot p — iuJkt— iat— T r t 



= ^ J d£ J kdk e 



Vl — — a — iT + iT r 

^iu> t— Tot— (r r — icr)e~ lkr t p i<^k* 



x 1 -p . / ~r --p \ _ ik . r (60) 

0J o - + li + [a + li r je lk r 

\2 p\ r. p—iujQt—TQt ^—iuj^t—icrt—Trt 



47r 2 7-i 7 



dj a — co»k — cr — ir + iT r 

giuj t— Tot— (F r — i<r)e~ lk ' r t pi^k^ 



x ; =^ — =— . (61) 

u - u k + «r + (or + «r r ) e - ik - r 



An approximate evaluation of the ^-integral followed by the further approximation ~ 
iir5(x), gives an estimation of the leading order contribution 

f i 1 (t)=T Q \fi(jt)+iu(t)] (62) 

fjL2(t) = TMt)+iv(t)] (63) 

where /i + zi/ is the complex-valued function 

u(t) + iu(t) =~ ^ : e -<« t-rot 

x(e -rot _ e -r r t^ e -r t-2^[r r -ia]t _ ^ty (64) 
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F. The master equation 



n T (t) 
PA L ) 


— e Pi\v) 




PoK 1 ) 


_ -2iuj o t-2T t I / n \ 

- e Po( u ) 




P-(t) 


— it,) f—'ZVnf-l-irrt-l-T f T /r^\ 




Piit) 


_ g — icuotSVot— iat—Trt pi 




PoU) 


= e -^*- rof +- +r "V5(o) + *(r + r r )!/(t) P i(o) + (r 


-r r )//(t) P i(o) 


P + o{t) 


= e --°*- r o*--- r -V5 (o) + (r + r r )/x(t)pi(o) + *(r 


-r>(Op 7 _(o) 


pXH) 


= e- 2r °*- 2r ^ P :(o) + (r + r r )«(<)^(0) 




p--(t) 


= e- 2r °*- 2r ^ P :(o) + (r - r P )«(t)^(o) 




p + -(t) 


= e 2iCT *- 2r °V^(0). 





Given the explicit form of the functions computed in Sec. II. E, we write the evolution 
equations in the following form 

(65) 
(66) 
(67) 
(68) 
(69) 
(70) 
(71) 
(72) 
(73) 

In the equations above we wrote the density matrix in a basis defined by \I), \0), |+), |— ), 
where 

|+> = ^(|01> + |10» (74) 
|-) = ^(|01)-|10». (75) 

We note that in this approximation, the diagonal elements of the density matrix propa- 
gator (i.e. the ones that map p%(0) to p&(£) where a,b G {O, +, — , J}) decay exponentially 
(except for the Pq which is functionally dependent on the other matrix elements due to 
normalization). We shall see in Sec. V that this behavior is in accordance with the Born- 
Markov approximation. The situation is different for the off-diagonal elements of the density 
matrix propagator, i.e. the ones that map p%(0) to Py(t) for a ^ a! and b ^ b'. They are 
given by more complex functions of time and they differ from the ones predicted by the 
Born-Markov approximation. 

As we have the solution p(t) = M t [p(0)], where M t is the density matrix propagator 
defined by Eqs. Q65H73p . we can identify the master equation through the relation pit) = 
MtiMf^pit)]]. We obtain 

Pi = -^op 1 ! (76) 
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pi = (-iu - 3r Q + iff + T r )p{ (77) 

p+ = (-iu - 3T - iff - r r )p+ (78) 

p^ = (-2iu o -2T )p I o (79) 

p = = (-iw - r + iff + r r ) Po + (r„ + r r )ai(t)pi + (r - r r )« 2 (t)pi (so) 

P + = = (-iw - r + iff + r r ) P + + (r + r r )« 3 (t)p^ + (r - r r )a A (t) P L (8i) 

p+ = -2(r + r r )p+ + (r + r> 5 (t)p{ (82) 

pz = -2(r - r r )p+ + (r - r r )« 6 (t)p{ (83) 

pi = 2(«T-r )^. (84) 



Explicit expressions for the functions of time atj(t), i = 1, . . . , 6 appearing in Eqs. (jTBHSlj) 
are given in the Appendix B. 

We see that the evolution equation for the reduced density matrix of the two qubits, 
while it is local-in-time, it does not have constant-in-time coefficients. Hence, it does not 
correspond to a Markov master equation of the Lindblad type. Again, we note that the 
non-Markovian behavior is solely found in the off-diagonal terms of the evolution law and 
that the diagonal ones involve constant coefficients. 

To facilitate comparison with the expressions obtained from the Born-Markov approxi- 
mations we cast equations (I76H841) into an operator form. 

'p = -i[H +Hi,p]+ ]T T lJ (S?S U) p + pS?S u) -2S^pS ( f > ) 

+ (T + r r )F t [p] + (r - r r )G t [p]. (85) 

The first term on the right-hand-side of Eq. (j85p corresponds to the usual Hamiltonian 
evolution: the total Hamiltonian is a sum of the free Hamiltonian Hq and of a dipole 
interaction Hamiltonian 

Hi = -ff(S- ®S + + S+®S-). (86) 
The second term in the right-hand-side of Eq. (I85p is the usual Lindblad term (see e.g., 



33j), where we defined T u = T 2 2 = T and T 12 = r 21 = T r . 

The last two terms contain effects that pertain to the off-diagonal terms of the reduced 
density matrix propagator and they are non-Markovian: F t and Gt are trace-preserving 
linear operators on the space of density matrices and their explicit form is given in Appendix 
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B. Assuming that the Markovian regime corresponds to the vanishing of F t and Gt, we find 
from Eqs. flB12f[BT3|) that in this regime the functions oti{t) in Eqs. (17614841) should reduce 



to the following constants 

oil — a i — 0; 0:3 = 05 = a§ = 2; a 2 = —2.. (87) 

In Sec. V, we shall discuss in more detail the physical origin of non-Markovian behavior 
in this two-qubit system. 



III. DISENTANGLEMENT OF TWO QUBITS 

In this section, we employ the results obtained above to study the evolution of the two 
qubits initially in an entangled state. We shall focus on the process of disentanglement 
induced by their interaction with the field. 



A. Class A states: Initial superposition of 1 00) and |11) 

We first examine the class of initial states we call Class A of the following type 



|Vo> = V 1 -P|00) + Vplll), (88) 

where < p < 1. Recall our definition |J) = |11) and \0) = 1 00) . From Eqs. (1281 - 1361) we 
obtain 

p{t) =p 2 e- ATot \I){I\ +e~ 2r ° t ^/p(l-p)(e 2iuJ ° t \I}(O\ + e- 2luJ0t \O}{I\) 
+Mt) - K 2 (t)}\-)(-\ + [«!(*) + Ka(t)]|+)<+| + [1 -p 2 e~ 4Tot - 2K 1 (t)]\0)(0\, (89) 



where the functions Ki(t) and K 2 (t) are given by Eqs. 

These results are quite different from those reported in Ref. which were obtained 

under the Born-Markov approximation. While the |/)(/| and \I)(0\ terms are essentially 
the same, the |— )( — | and |+)(+| ones are not, as they involve non-diagonal elements of the 
density matrix propagator. For comparison, we reproduce here the explicit form of these 
matrix elements in our calculation 

pX = Pf^r; e ~ 2rot ^ rot - e ~ Trt ? ^ 

p - = pe ~2r t ( l-r t_ e ~rrty^ (91) 
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FIG. 1: A plot of the concurrence as a function of time for initial Class A state (I88p with p = 0.5 
for different values of uj Q r. 



and in that of Ref. |40j (translated into our notation): 

Pt(t) 



T — T r 

p-jt) = p £^zIl e -2rot fe 2r ri 



-2F t 



)■ 



(92) 
(93) 



r + r r 

For large values of r, T r « T and the expressions above coincide. However, for oj Q r « 1 
their difference is substantial. In this regime, To — r r and at times Tot ~ 1 we obtain 
(r — T r )t << 1. According to the Markovian results of 40j, in this regime the |+)(+| 
term is of order O(A ) and hence comparable in size to the other terms appearing in the 
evolution of the density matrix. However, according to our results, which are based on the 
full non-Markovian dynamics, the |+)(+| term is of order r °^ r and hence much smaller. 

In general, for uor << 1, we find that the |— )(— | and |+)(+| terms contribute little to 
the evolution of the reduced density matrix and they can be ignored. Since these terms are 
responsible for the sudden death and subsequent revival of entanglement studied in [40J], we 
conclude that these effects are absent in this regime. Indeed, this can be verified by the study 
of concurrence as a function of time as appearing in Figs. 1 and 2. For large inter-qubit 
separations and for specific initial states {p > 0.5) there is sudden death of entanglement-but 
no subsequent revivals. 
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FIG. 2: A plot of the concurrence as a function of time for initial Class A state (|88p for different 
values of p and fixed inter-qubit distance (oj r = 1). 

B. Initial antisymmetric |— ) Bell state 



We next consider the case of an initial antisymmetric |— ) state for the two qubits. We 
find 

p(t) = e-^'-^HH + (l -e- 2 ^- 1 ^) \0)(0\ (94) 

We see that the decay rate is r — T r . The effect of photon exchange between the qubits 
essentially slows down the overall emission of photons from the two-qubit system (sub- 
radiant behavior). As the qubits are brought closer, the decay is slower. At the limit r — ► 0, 
there is no decay. 

n 

The results agree qualitatively with those obtained in Ref. [33j through the Born-Markov 
approximation. Fig. 3 shows the time evolution of concurrence for an initial antisymmetric 
|— ) state and for different inter-qubit distances. 



C. Initial symmetric |+) Bell state 



For an initial symmetric |+) Bell state the reduced density matrix of the two qubits is 



p(t) = e-^+^l+X+l + (l - e - 2 ^ +T ^) \0)(0\. 



(95) 
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FIG. 3: A plot of the concurrence as a function of time for the initial antisymmetric |— ) Bell state 
and for different values of the inter-qubit distance r (in units of The decay of concurrence 

proceeds more slowly when the qubits are closer together. 




r t 



FIG. 4: Plot of the concurrence as a function of time for initial symmetric |+) Bell state and for 
different values of the inter-qubit distance r (in units of lo^ 1 ). 



Here we have a super-radiant decay with rate given by T + r r . Note that for both the 
antisymmetric |— ) and symmetric |+) states, a resonance appears when r = mru^ 1 , and the 
decay becomes radiant. Fig. 4 shows the behavior of concurrence, in qualitative agreement 
with the corresponding results of Ref. j^. 
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IV. DECOHERENCE OF QUBITS 



In this section we study the evolution of factorized initial states, in order to understand 
how the decoherence of one qubit is affected by the presence of another. 

A. One qubit in vacuum state 

We assume that the first qubit is prepared initially in a superposition of the |0) and |1) 
states, and that the second qubit lies on the ground state. Therefore, the initial state is of 
the form 



y/p\l) +yjl- P \0)J ® |0> = y/p\10) +y/l- p|00>. (96) 

From Eqs. ( I2"8"H3"6"1) . we obtain the density matrix of the combined qubit system 

p{t) =p(|w + | 2 |10)(10| + |v_| 2 |01)(01| + u*v+|01)(10| + u_i£|10)(01|) 
+ y/p{l-p) (<|10)(00| +«1|01)<00| + u+|00)(10| + v_|00)(01|) 

+ (l-p(K| 2 + M 2 ))|00>(00|, (97) 



where the functions v±(t) are given by Eq. (jHJ). 

The two qubits become entangled through their interaction via the EM field bath. To 
study the decoherence in the first qubit, we trace out the degrees of freedom of the second 
one, thus constructing the reduced density matrix p 1 



h{t) = PK| 2 |1>(1| + yjp(l-p(v+\0)(l\ + <|1)0|) + (l - P \v + \ 2 ) |0)(0|. (98) 

At the limit of large interqubit distances, T r = = c(r), whence v + ~ e -iu o t-v t^ 
off-diagonal elements decay within a characteristic time-scale of order IV . These results 
coincide with those for the single qubit case-see Refs. 34, 3s| . However, in the regime 
uj r « 1, the results are substantially different. The entanglement with the second qubit 
leads to a departure from pure exponential decay. In this regime, T r ~ IV This implies 
for times longer than IV 1 a substantial fraction of the off-diagonal elements persists. This 
decays eventually to zero within a time-scale of order [To — r r ] _1 >> IV ■ Hence, the qubit 
preserves its coherence longer. (At the limit r — > there is no decoherence.) 

The reduced density matrix of the second qubit is 



k(t) = pM 2 |i>UI + v^(i-p)(*>-|o>(i| + u_|l>0|) + (l ~p\v-\ 2 ) |0)(0|. (99) 
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Note that at the limit of small inter-qubit distances the asymptotic behavior ( for Tot » 
1) of Pi(t) is identical to that of p 2 (t). The second qubit (even though initially on its ground 
state) develops a persistent quantum coherence as a result of the interaction with the first 
one. 



B. One qubit in excited state 

We also consider the case of a factorized initial state, in which the second qubit is excited 



y/p\l) + yjl -p\0)j ® |1> = y/p\I) + yj\ -p|01). (100) 

This system behaves differently from that of Sec. IV. A. The matrix elements of the reduced 
density matrix read 

p\l = pe-^ (101) 



Pl\ = ^p(l- P )e- 2 ^ t - 2Tot v*_ (102) 
PoS = Jp(l-p)Pi(t) (103) 



pl° = y/p(l-p)Mt) (104) 

Pol = (l-p)|«+| 2 +P«i (105) 
p% = {l-p)v + v*_+ P K 2 , (106) 

where the functions /ii )2 are given by Eqs. (|62H64|) . the functions by Eqs. (!56|l and 
f l57|) and the functions v± by Eq. (jUj) . 

The reduced density matrix of the first qubit is 

p x (t) = ((1 -p)\v + \ 2 + P k, +pe- 4Tot ) |1)(1| + (l - (1 -p)\v + \ 2 - V k x -pe- 4r °') |0)(0| 

+ y/p(l-p)((p* 2 + e- 2i ^- 2rot <)|0)( 1 l + (P2 + e 2 ^*- 2rot <)|l}(0|)(107) 

At times t » Tq 1 , the decay of the off-diagonal elements at is dominated by the function 
p 2 , which then reads 

"»« - r. + ^w PS)''"" " e " 0<) ' (1 ° 8) 

In the limit of large inter-qubit distance r, the off-diagonal element falls like r _2r °*, while 
in the small r limit it falls like e~ r °*. Comparing with the case of Sec. IV.A, we see that the 
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initial excitation of the second qubit leads to a lesser degree of entanglement of the total 
system, as it cannot absorb any virtual photons emitted from the first one. For this reason, 
the asymptotic decoherence rate does not vary significantly with the interqubit distance. 



V. DISCUSSION 



In this section, we first summarize our results for the evolution of various initial states, we 
then discuss the origin of the non-Markovian behavior and finally, a possible limitation of our 
results that is due to the restricted domain of validity of the Rotating- Wave approximation. 



A. Description of the results 

For an initial Class A state (1551) . the \I) component decays to the vacuum, but it also 
evolves into a linear combination of antisymmetric |— ) and symmetric |+) Bell states. How- 
ever, if the qubits are close together the evolution to Bell states is suppressed. This behavior 



is qualitatively different from that of Ref. [40|, which was obtained through the Born-Markov 
approximation. The corresponding terms differ by many orders of magnitude at the physi- 
cally relevant time-scales. As a consequence, we find that there is neither sudden death nor 
revival of entanglement in this regime. 

In retrospect, this difference should not be considered surprising. The Born-Markov 
method involves two approximations: i) that the back-action of the field on the atoms is 
negligible and ii) that all memory effects in the system are insignificant. When the qubits 
are found within a distance much smaller than the one corresponding to their characteristic 
wavelength, it is not possible to ignore back-action. The virtual photons exchanged by the 
qubits (at a rate given by r r ) substantially alter the state of the electromagnetic field. 

On the other hand, the effect of virtual photons exchange between qubits drops off quickly 
at large separations r - the two qubits decay almost independently one of the other. Hence, 
the Born-Markov approximation - reliable for the case of two separate qubits each interacting 
with an individual field - also gives reasonable results for the two qubits interacting with a 
common field. In this regime sudden death is possible, but not revival of entanglement. In 
this sense, our results effectively reduces to those of Ref. 



: when the distance between 



the qubits is much larger than any characteristic correlation length scale of the system it 
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looks as though the two qubits are found in different reservoirs. 

Therefore, our results for initial states of Bell type are qualitatively similar to those 
obtained in Ref. 33| through the Born-Markov approximation. The symmetric |+) state 



decays super-radiantly and the antisymmetric |— ) one sub-radiantly. For small values of 
u r, the antisymmetric |— ) state decays very slowly and entanglement is preserved. 

Concerning decoherence, when the inter-qubit separation is large and the second qubit 
lies on its ground state, our two qubit calculation reproduces previous results on the single 



qubit case [34|, |35j. However, if the qubits are close together the coherence is preserved 



longer. These results seem to suggest that in a many-qubit system, the inter-qubit quantum 
coherence can be sustained for times larger than the decoherence time of the single qubit 
case. This may suggest some physical mechanisms to resist decoherence in multi-qubit 
realistic systems. 



B. The origin and significance of the non-Markovian behavior 

With the exact solutions we have obtained (under weak coupling but no Born or Markov 
approximation) for the two qubit - EMF system we want to elaborate on the origin and 
significance of the non-Markovian behavior started in the Introduction. In the evolution 
equations ( I65T473I) we note that the diagonal terms of the reduced density matrix propagator 
all decay exponentially. Their decay rate is therefore constant. It is well known that this 
feature is a sign of Markovian behavior. In fact, it characterizes the domain of validity 
of Fermi's golden rule: one could obtain the decay rates by direct application of this rule. 
Hence, as far as this part of the evolution is concerned, our results are fully compatible with 
the Markovian predictions. 

However, the behavior of the non-diagonal terms in the reduced density matrix propagator 
is different. A look at Eqs. ( 16511731) will convince the reader that the only non-zero such 
terms are ones that describe the effect of successive decays, for example the |11) state first 
decaying into |— ) and then |— ) decaying into the ground state |00). Hence, the pZ{t) term 
consists of one component that contains the remaining of the |— )(— | part of the initial state 
and another component that incorporates the decay of the |11)(11| part of the initial state 
towards the state |— ). In our calculation, the latter term is encoded into the functions ni^{t), 
which are obtained by squaring the amplitudes Vk(t) as in Eqs. (139H401) . At second loop 
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order, these amplitudes are obtained from the summation of two Feynman diagrams -see Eq. 
(|5ip . The structure of the poles in Eq. floTj) reveals that the first Feynman diagram describes 
the decays of the |11) state, while the second one corresponds to processes involving the 1 01) 
and 1 10) states. When we construct the evolution functions ^1,2(^)5 we obtain terms that are 
both diagonal and off-diagonal with respect to the two types of process. It is the presence of 
the off-diagonal terms that is primarily (but not solely) responsible for the deviation of our 
results from the Markovian prediction. In the Markov approximation, the corresponding 
term involves summation (subtraction) of probabilities rather than amplitudes. 

To justify the last statement, we note from Eqs. ( I28H361) that the evolution decouples 
the diagonal from the off-diagonal elements of the density matrix (the Markovian equations 
also have this property). Hence, the probabilities p a (t) = p^(t),a G {/, +, — ,0} evolve 
autonomously. Time-homogeneity (i.e. Lindblad time evolution) implies that their evolution 
can be given by a transfer matrix 

Pa(t) = Y,T b aPb (t). (109) 

b 

Here T fc a is the constant decay rate for the process b — > a. Noting that Eqs. (12814361) contain 
no transitions between + and — and no transitions from — to J, Eq. (I109p yields 

p_(t) = -2(r -r r ) P _(t) + w Pl (t), (110) 

where w = Tf . Since pi(t) is determined by Eq. (ITS"]) as p/(0)e~ 4r °*, we obtain for the initial 
condition p_(0) = 0, 

p_(f) ~ j9 / (0)(e- 2(r °- r " ) * - e" 4rot ), (111) 

in full agreement with Eq. ( 193]) obtained from the Lindblad master equation. 

The derivation of Eq. (II lip provides an example of a more general fact: the Markovian 
assumption forces the off-diagonal terms of the density matrix propagator (in this case the 
one mapping pj(0) to pZ{t)) to be subsumed by the diagonal ones. The off-diagonal ele- 
ments can have no independent dynamics of their own, unlike what would be the case if 
they were derived from a full calculation. One should also note that the diagonal terms 
of the propagator are obtained from the first order perturbation theory. Since they deter- 
mine the off-diagonal terms within the Markov approximation, the latter only contains the 
information obtained from first-order perturbation theory. However, in our calculation the 
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off-diagonal terms involve second-order effects and hence, reveals dynamical correlations 
that are inaccessible within the context of the Born-Markov approximation. 

For the reasons above, the matrix elements (JH2HS2D given by 40) obtained under the 
Markovian approximation-carry the characteristic superradiant and subradiant behavior of 
the decay of the |+) and |— ) states respectively, a property that does not arise from our 
calculation. 

Moreover, we note that Eq. (IllOp obtained by making the Markov approximation is 
related to arguments from classical probability, i.e. it involves addition of the probabilities 
associated to different processes 52|. On the other hand, a proper quantum mechanical 
calculation involves the addition of amplitudes-e.g. the ones appearing in the definition of 
the functions and as such it must also contain 'interference' terms. In effect, the 

Markovian approximation introduces by hand a degree of partial 'decoherence' (or classical- 
ity) and we think that this explains why in general it predicts a faster classicalization of the 
system than the full analysis does. 

We believe that the feature discussed above is generic. In the single qubit system, it was 
not present because its structure was too simple. There could be no intermediate decays. 
However, this effect should in principle be present in any system that contains intermediate 
states. The Markovian approximation would then be valid only if specific conditions hold 
that render the 'interference' terms negligible-for example if there exists a sharp separation 
of the relevant timescales. 

To summarize, the Markov approximation essentially ignores interference terms that are 
relevant to processes that involve successive decays. These processes appear through off- 
diagonal terms in the reduced density matrix propagator. The Markov approximation mis- 
represents the intrinsic dynamics for these terms and ties them -by forcing additivity of 
probabilities-to the evolution of the diagonal ones. As a results, the off-diagonal terms are 
subsumed by the diagonal ones. 



C. The use of the rotating wave approximation 

Finally, we add a few words on the accuracy of our model for the interaction between 
the 2LA's and the EMF. The Hamiltonian (j2J) is obtained from the study of the interaction 
of the atomic degrees of freedom with the electromagnetic field. Its derivation involves 
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the dipole and the rotating wave (RW) approximations-see Appendix A in Ref. 34|. The 
RW approximation consists in ignoring rapidly oscillating terms in the interaction-picture 
Hamiltonian and keeping only the ones that correspond to resonant coupling. (One ignores 
processes during which a photon is emitted and the atom becomes excited.) The terms that 
are dropped out in the RW approximation oscillate with a frequency of order uo . 

For a single qubit system the RW approximation is self-consistent. However, in the 
two-qubit system we keep terms in the Hamiltonian that vary in space as e* kr . For the RW 
approximation to be consistent, one has to assume that k • r << uj t. Since k is peaked 
around the resonance frequency, this condition is equivalent to r << t. Since the physically 
interesting time-scale for the study of disentanglement and decoherence corresponds to 
Tq 1 ~ [A 2 co> ] -1 , we expect the RW approximation to be reliable in this context, as long 
as \ 2 u r « 1. This is sufficient for realistic situations, in which r is bounded by the 
size of the cavity and A 2 << 1. However, the condition above does not hold at the 
formal limit r — > oo. In this regime the RW approximation may break down. Indeed, in 
section IV. B the reduced density matrix for the single qubit in the limit r — > oo does not 
coincide with the corresponding expression obtained in the study of the single qubit. The 
presence of an excited second qubit, even if it is situated far away, affects the time evolution 
significantly. This effect is also present in the Born-Markov approximation. This is arguably 
an unphysical behavior, and we believe that it arises as an artefact of the RW approximation. 
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APPENDIX A: PERTURB ATIVE RESUMMATION 



We notice that the action of the operator [(E — Ho)~ l Hi\ 2 on any linear combination 
of the vectors |0; 0, 1) and |0; 1, 0) yields another linear combination of these vectors. If 
we denote by p n and q n the coefficients of \0; 0, 1) and \0; 1, 0) respectively after the n-th 
action of [(E — iJo) 1 -^/] 2 (2n-th order of perturbation theory) we obtain 

/ 




1 



E-n 



a(E) (3*{E,r) 
P(E,r) a(E) 



--: A 



\ Qn-l J 



\ Qn-l J 



(Al) 



Consequently, 





= A n 


("") 


\q n ) 







(A2) 



and summing over all n we obtain 




n=0 





(A3) 



To compute (z; 0, 1\(E - //) _1 |0; 01) and (z; 1, 0\(E - //) _1 |0; 01) according to the per- 
turbation series we have po = (E — c^ ) _1 , qo = 0. The results in the main text then follow. 

The resummation for (z; 0, 0\(E — if ) _1 |0; 0, 1) proceeds from the fact that the terms in 
the 2n + 1 order of the perturbative expansion are of the form 

g k ff e ik-r/2 



-b{\0, 0,0) 



(A4) 



We then note that 



_ l a(E)+(3(E,r)e i ^ ^ n _ 1 

— [ J? ~ IS] 



, , k ■ (A5) 

E — to o 

where Ck — (E — ^ ) _1 - This geometric series can be summed to give the result quoted in 
the main text. 

A similar summation is achieved for the terms (z; 0, 1\(E — H)~ l \0; 1, 1) and (z; 1, 0\(E — 
iJ) 1 10; 1, 1). In the 2n-th order of perturbation theory the action of the series on O; 1, 1) 



yields 



E 



9k 



(E -2fl)(E -fi-Wk) 



(e- ik ' r/2 7™|0; 0, 1) + e ik ' r/2 ^|0; 1, 



(A6) 
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APPENDIX B: THE TIME-DEPENDENT TERMS OF THE MASTER EQUA- 
TION 



The time- dependent coefficients a iy i = 1, . . . 6 appearing in the master equation 
are defined as 



ai(t) 
a 2 (t) 
a 3 (t) 

a 5 (t) 
a 6 (t) 



Using the notation 



d 



(M*) 



t+r t-icrt~T r t \ e 2T t+2icrt+2r r t 



e iiOot+T Q t-iat-r r t\ e 2T t 



dt 

d_ 

dt 

— (n(t)e iu ' ot+Tot+iat+T ' rt ^j e 2r °* 
dt 

ij t (v(t)e iuJ ° 



t+r t+iat+Trt) e 2T t-2i(rt-2V r t 



t+2T r t\ e 2V Q t-2V r t 



A +I = 




Aoi = 


\0)(I\ 


A ++ = 


\+)(+\ 


A+- = 


\+)(-\ 


= 


l-X-l 



the linear functional F t and G t of Eq. flSBT) are defined as 



(Bl) 
(B2) 
(B3) 
(B4) 
(B5) 
(B6) 



(B7) 
(B8) 
(B9) 
(BIO) 
(Bll) 



F t [p] = -[2 - a 5 (t)}(A+ipA-ii - AoipA j OI ) - [2 - a 3 (t)]A OI pA ++ - [2 - al(t)]A ++ OI 

+ a 1 (t)i jM+- + a 1 (£)it-Mo/ ( B12 ) 

G t [p] = -[2 - a 6 (t)](i_ / pi t _ / - A OI pA ] OI ) + [2 + a 2 {t)]A OI pA^ + [2 + a* 2 (t)]A^0 OI 

+ a 4 (t)AoipA j + _ + al{t)A ] + _ P A ] OI . (B13) 



APPENDIX C: A GRASSMANN PATH-INTEGRAL DERIVATION 



In Refs. 34], |35| the evolution of the reduced density matrix for a 2LA interacting 
with the EM field was determined with a version of the Feynman- Vernon path integral 
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method, using Grassmann variables for the atomic degrees of freedom. The same method 
can be applied to the present problem, even though the perturbative approach turned out 
to be more convenient because the vacuum initial state allowed for a summation of the 
perturbative series. However, the path-integral method may allow for a simpler treatment 
of other systems, such as an N-qubit system or the EM field at a finite temperature. For 
this reason, we include a sketch of the path-integral treatment in this Appendix, noting that 
it leads to the same results as the ones obtained in Sec. II. 



1. Coherent state path integral 

The coherent states for the EMF (bosonic) and 2LA (qubit) states are defined respectively 

by 

\z k ) = exp(^ k 6 k )|0 k ) (CI) 
\rj^) = expfasi B) )|0) (B) . (C2) 

The states |0 k ) and |0)( n ) are ground states of the electromagnetic field's k^ 1 mode in 

th 

vacuum and the n qubit in its lower state, repectively. The bosonic coherent states are 
labeled by a complex number, z k , and the qubit coherent states are labelled by an anti- 
commuting (i.e. Grassmannian) number, rf n \ A non-interacting eigenstate of two qubits 
and the electromagnetic spectrum can be written as the direct product of coherent states, 

\{z k },V W ,V {2) ) = IK}) ® \v m ) ® \v {2) )- (C3) 

In a path integral approach the transition amplitude between some chosen initial and 
final states is divided into many infinitesimal time steps. The resulting path integral corre- 
sponds to the matrix element of the evolution operator. In the coherent-state path integral 
representation the Hamiltonian is given by, 

exp^yw + ffarfW + £ k z k z' k ] nUJ ° r,rj +n \ 



+ E (AL n V n) 4 + A k "W W )(C4) 



Spatial dependence is carried in the coupling constants, 



\(n) 

\W = ^ e ~^ (C5) 

C<J k 
\(n) 

- X (n) = ^_ e +ik.r n (C6) 
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» 



Pow = 



/i+l,kl — (1 - ^kC)/jkl - «e Em ^'k^ S'jql' /okl — <^kl- 



The transition amplitude for N atoms after the j-th infinitesimal time steps is, 
K(je, 0) = exp { £ [ £ + £ $ J + £ \ £ /i^oi + E } • 

^ n L m kl k 1 fw 

(C7) 

The functions in the path integral are determined by finite-difference equations 

= (1 - iu e)^f m) - ie E k 4 nm) = 4 



(C8) 



(C9) 



2. Reduced Density Matrix 

The reduced density matrix of the two 2LAs is 

pit) = J dn(rj?> )d//(rf 2) )^tf ) )d / i(rf 1) )p(0)^(t,0). (CIO) 
The dynamics of this open system is described by the density matrix propagator 

Mt, o) | T=0 = exp { e k n) m (nm) vt ] + nT^'{t) {nm) 4 m) + 4 n) E mMU m) 

(Gil) 

A general initial two-qubit density matrix written in the Grassmann representation is: 

f pB8 p8? p?S pSS \ ^ i ^ 



p(o) = (i # # 



Pffi P°o\ P?o Pfi 



Poo Poi PlO Pll 



V Poo Pol p\1 p\l J \ % ( No (2) ) 



1(2) 



(C12) 



In addition to Eqs. Eq. (1C8IIC9l) . we derive finite difference equations for the expanded 
functions appearing in the reduced density matrix elements 



CiCi = a - 2^ e)< b vr m} Agvsvr - * e s w k vr j ( ci3 ) 



^j+^+i.k = (1 - iu; e - ^k)V>. 



(n) , (nm) , (aft) 



£^ 



(C14) 



= (! - ^ e - - <c E AfkVrVi? - ie E ^Vf V? (C15) 
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For a two-atom system, it is convenient to represent the functions appearing in the propa- 
gator by column vectors, 



^,(22) \ 

^(21) 
^(12) 

^(22) \ 



22) 



^(21)^(12) 



(C16) 
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( 4M ] 
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At the continuous limit the finite difference equations yield a set of differential equations for 
the column vectors above. 
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\A™ 

The solution of Eqs. (1C18I4U221 allows one to fully reconstruct the propagating amplitude 
and from this the elements of the reduced density matrix for the qubits. These equations 
can be solved by implementing an approximation scheme similar to that employed for the 
operator method in the main text. We do not provide the details, but only note that the 
results of the two methods coincide. 
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